library(here)
#analisis
library(ggsignif)
library(ggrepel)
library(ggpubr)
library(inlmisc)
library(nortest) #para testear distribucion
library(skimr) #provides a frictionless approach to summary statistics w
library(easystats) # multiples unciones analiticas
library(lme4)
library(skimr)
library(readxl)
# vizualizacion
library(ggridges)
library(sf)
library(GGally)
library(tidyverse, quietly = TRUE)
library(knitr, quietly = TRUE)
library(kableExtra)
library(raster)
library(egg)
library(car) #Variance inflation Factor
library(ggthemes)
library(sjPlot)Análisis espacio-temporal datos histórica artesanal de la pesquería de Bacalao
Analisis exploratorio 1986-2021
ANTECEDENTES
El presente reporte tiene un codigo autocontenido con los Analísis Exploratorios de Datos (AED) de la pesquería de bacalao de profundidad extraído por la flota pesquera artesanal (espinel). Los analisis aca descritos estan compuestos de dos bases. La primera es de los registros oficiales de captura recogidos por el SERNAPESCA. En segundo lugar, se presentan el AED de los datos del monitoreo d la pesqería llevado a cabo por el IFOP a través del Programa de Seguimiento de Pesquerías Demersales del Departamento de Evaluación de Pesquerías.
El principal objetivo es identificar vacios y fortalezas de los datos en terminos espacio-temporales, y con ello tomar decisiones para el paso posterior de modelación de la dinámica poblacional del recurso para realizar asesoría en terminos de estatus y recomendación de una Captura Biologicamente Aceptable (CBA).
ANÁLISIS EXPLORATORIO DE DATOS (AED)
Cargo librerías necesarias para el análisis exploratorio de los datos de las distintas bases de bitácora, tallas y biológico.
# A function for dotplots
multi_dotplot <- function(filename, Xvar, Yvar){
filename %>%
ggplot(aes(x = {{Xvar}})) +
geom_point(aes(y = {{Yvar}}),
alpha=0.4) +
theme_bw() +
coord_flip() +
labs(x = "Order of Data")}Identifico los directorio de trabajo y leo las tres bases de datos; a saber:
- Datos FIP (Rendimiento) 96-32
- Bitácora
- Biologico
- Tallas
- Desembarques (1985-2022)
bit <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BITACORAS ESPINEL bacalao_corr.xlsx", sheet = "BITACORAS_ESPINEL_bacalao", guess_max = 100000)
long <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/LONGITUD_ESPINEL_bacalao.xlsx")
bio <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BIOLOGICO ESPINEL bacalao.xlsx",
sheet = "BIOLOGICO_ESPINEL_bacalao")
landing <- read_excel("Artesanal_Historica/DESEMBARQUE HISTÓRICO.xlsx",
skip = 2)
bit8696<- read_csv2("bit_art_86_20.csv")Identificamos los registros asociados recurso objetivo
- Codigo del Recurso: 37
- Nombre común: Bacalao de profundidad
- Nombre científico: Dissostichus eleginoides
Idetifico los registros por recurso y por base para luego filtrar.
dim(bit)
[1] 14208 77
table(bit$COD_ESPECIE)
2 4 5 6 14 15 16 17 18 20 22 23 24
7 1 14 25 147 6 89 16 698 1 25 3 7
25 27 33 35 37 42 43 49 53 56 59 74 75
3 11 2 2 11074 3 9 21 1 5 35 1 28
81 88 99 103 104 105 106 108 119 127 129 144 194
19 1 159 21 24 2 7 44 1 1 1 1 91
199 212 226 239 243 323 324 325 334 844 849 920 984
178 18 1 22 1 2 50 6 830 3 138 1 1
989 999 1309
36 249 58
dim(long)
[1] 15831 24
table(long$COD_ESPECIE)
37
15831
dim(bio)
[1] 127409 36
table(bio$COD_ESPECIE)
37
127409
# las bases de biologicos y longitud solo teienen registros de bacalao.Comenzar a trabajar bases por separado
DESEMBARQUES
Tabla con los desembarques oficiales (Sernapesca, 2022)
kbl(landing, booktabs = T,format = "html",
caption = "Desembarque Bacalao Artesanal por Región") %>%
kable_styling(latex_options = c("striped",
"condensed","scale_down"),
full_width = FALSE) | ...1 | AyP | TPCA | ANTOF | ATCMA | COQ | VALPO | LGBO | MAULE | ÑUBLE | BBIO | ARAUC | RIOS | LAGOS | AYSEN | MAG | ...17 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1985 | NA | 8 | 102 | 146 | 162 | 1674 | NA | 349 | NA | 1527 | NA | NA | 63 | 28 | NA | 4059 |
| 1986 | NA | 842 | 1531 | 213 | 329 | 1242 | NA | 88 | NA | 1546 | NA | NA | 311 | 6 | NA | 6108 |
| 1987 | NA | 182 | 335 | 171 | 51 | 671 | NA | NA | NA | 1264 | NA | NA | 689 | 21 | NA | 3384 |
| 1988 | NA | 131 | 174 | 251 | 42 | 670 | NA | 46 | NA | 1597 | NA | NA | 877 | 8 | NA | 3796 |
| 1989 | NA | 217 | 175 | 177 | 109 | 756 | NA | 114 | NA | 1859 | 10 | NA | 1462 | 8 | NA | 4887 |
| 1990 | NA | 335 | 161 | 126 | 100 | 444 | NA | 5 | NA | 2620 | NA | NA | 1789 | 36 | NA | 5616 |
| 1991 | NA | 280 | 186 | 142 | 223 | 519 | NA | NA | NA | 1901 | NA | NA | 680 | NA | NA | 3931 |
| 1992 | NA | 30 | 45 | 77 | 53 | 477 | NA | 131 | NA | 2108 | 1 | NA | 741 | 1 | NA | 3664 |
| 1993 | NA | 37 | 4 | 66 | 81 | 1676 | NA | 295 | NA | 1111 | NA | NA | 807 | NA | 45 | 4122 |
| 1994 | NA | 24 | 139 | 213 | 68 | 671 | NA | 526 | NA | 1920 | NA | NA | 1374 | 41 | 411 | 5387 |
| 1995 | NA | 252 | 259 | 286 | 65 | 206 | NA | 190 | NA | 1347 | NA | NA | 1680 | 31 | 266 | 4582 |
| 1996 | NA | 476 | 220 | 121 | 177 | 272 | NA | 226 | NA | 1390 | 2 | NA | 1989 | 22 | 83 | 4978 |
| 1997 | NA | 36 | 82 | 47 | 57 | 98 | NA | 116 | NA | 1077 | 1 | NA | 1803 | 22 | 83 | 3422 |
| 1998 | NA | 9 | 67 | 169 | 61 | 155 | NA | 399 | NA | 1284 | NA | NA | 2038 | NA | 11 | 4193 |
| 1999 | NA | 42 | 72 | 168 | 87 | 484 | NA | 406 | NA | 1369 | NA | NA | 2910 | NA | 270 | 5808 |
| 2000 | NA | 285 | 98 | 153 | 105 | 200 | NA | 292 | NA | 800 | NA | NA | 3515 | NA | 345 | 5793 |
| 2001 | NA | 198 | 85 | 127 | 113 | 191 | NA | 136 | NA | 1049 | NA | NA | 2045 | NA | NA | 3944 |
| 2002 | NA | 154 | 44 | 97 | 57 | 188 | NA | 186 | NA | 732 | NA | NA | 3097 | NA | 10 | 4565 |
| 2003 | NA | 58 | 25 | 51 | 27 | 338 | NA | 123 | NA | 752 | NA | NA | 3189 | NA | 179 | 4742 |
| 2004 | NA | 76 | 32 | 63 | 23 | 308 | NA | 132 | NA | 412 | NA | NA | 2303 | NA | 70 | 3419 |
| 2005 | NA | 111 | 14 | 47 | 26 | 84 | NA | 175 | NA | 475 | NA | NA | 2021 | NA | 325 | 3278 |
| 2006 | NA | 87 | 3 | 50 | 50 | 92 | NA | 159 | NA | 243 | NA | NA | 1367 | NA | 40 | 2091 |
| 2007 | 19 | 43 | 7 | 22 | 25 | 43 | NA | 84 | NA | 189 | NA | 501 | 1157 | NA | NA | 2090 |
| 2008 | 4 | 48 | 21 | 9 | 3 | 32 | NA | 70 | NA | 219 | NA | 349 | 803 | NA | NA | 1558 |
| 2009 | 13 | 85 | 24 | 20 | 16 | 79 | NA | 112 | NA | 239 | NA | 522 | 571 | NA | NA | 1681 |
| 2010 | 13 | 43 | 13 | 58 | 11 | 37 | NA | 63 | NA | 109 | NA | 465 | 655 | NA | NA | 1467 |
| 2011 | 19 | 82 | 20 | 46 | 6 | 49 | NA | 90 | NA | 276 | NA | 613 | 987 | NA | NA | 2188 |
| 2012 | 25 | 68 | 11 | 82 | 14 | 26 | NA | 105 | NA | 317 | NA | 290 | 1126 | NA | NA | 2064 |
| 2013 | 26 | 58 | 12 | 68 | 23 | 52 | NA | 214 | NA | 186 | NA | 277 | 622 | NA | 20 | 1558 |
| 2014 | 13 | 48 | - | 59 | 34 | 64 | NA | 171 | NA | 129 | NA | 252 | 453 | NA | 57 | 1280 |
| 2015 | 6 | 52 | 13 | 53 | 38 | 93 | NA | 119 | NA | 261 | NA | 298 | 638 | NA | 38 | 1609 |
| 2016 | 17 | 51 | 20 | 59 | 16 | 70 | NA | 71 | NA | 335 | NA | 327 | 755 | NA | NA | 1721 |
| 2017 | 20 | 92 | 33 | 89 | 33 | 117 | NA | 116 | NA | 316 | NA | 312 | 790 | NA | 28 | 1946 |
| 2018 | 29 | 78 | 5 | 74 | 49 | 78 | NA | 113 | NA | 251 | NA | 218 | 620 | NA | 25 | 1540 |
| 2019 | 26 | 55 | 7 | 63 | 35 | 42 | NA | 102 | NA | 248 | NA | 203 | 696 | NA | 94 | 1571 |
| 2020 | 5 | 18 | 7 | NA | 1 | 27 | NA | 78 | NA | 145 | NA | 229 | 261 | NA | 32 | 803 |
| 2021 | 2 | 83 | 10 | 138 | 17 | 52 | NA | 206 | NA | 304 | NA | 500 | 517 | NA | 21 | 1850 |
| 2022 | 26 | 100 | - | 148 | 58 | 75 | - | 143 | - | 318 | - | 332 | 754 | - | 29 | 1983 |
Ploteo los desembarques por region y por año
landing <- as.data.frame(lapply(landing, as.double))
str(landing)
'data.frame': 38 obs. of 17 variables:
$ ...1 : num 1985 1986 1987 1988 1989 ...
$ AyP : num NA NA NA NA NA NA NA NA NA NA ...
$ TPCA : num 8 842 182 131 217 335 280 30 37 24 ...
$ ANTOF: num 102 1531 335 174 175 ...
$ ATCMA: num 146 213 171 251 177 126 142 77 66 213 ...
$ COQ : num 162 329 51 42 109 100 223 53 81 68 ...
$ VALPO: num 1674 1242 671 670 756 ...
$ LGBO : num NA NA NA NA NA NA NA NA NA NA ...
$ MAULE: num 349 88 NA 46 114 5 NA 131 295 526 ...
$ ÑUBLE: num NA NA NA NA NA NA NA NA NA NA ...
$ BBIO : num 1527 1546 1264 1597 1859 ...
$ ARAUC: num NA NA NA NA 10 NA NA 1 NA NA ...
$ RIOS : num NA NA NA NA NA NA NA NA NA NA ...
$ LAGOS: num 63 311 689 877 1462 ...
$ AYSEN: num 28 6 21 8 8 36 NA 1 NA 41 ...
$ MAG : num NA NA NA NA NA NA NA NA 45 411 ...
$ ...17: num 4059 6108 3384 3796 4887 ...
summary(landing)
...1 AyP TPCA ANTOF
Min. :1985 Min. : 2.00 Min. : 8.00 Min. : 3.00
1st Qu.:1994 1st Qu.:11.25 1st Qu.: 44.25 1st Qu.: 12.75
Median :2004 Median :18.00 Median : 77.00 Median : 32.50
Mean :2004 Mean :16.44 Mean :128.26 Mean : 112.67
3rd Qu.:2013 3rd Qu.:25.25 3rd Qu.:148.25 3rd Qu.: 111.25
Max. :2022 Max. :29.00 Max. :842.00 Max. :1531.00
NA's :22 NA's :2
ATCMA COQ VALPO LGBO
Min. : 9.0 Min. : 1.00 Min. : 26.0 Min. : NA
1st Qu.: 58.0 1st Qu.: 23.50 1st Qu.: 65.5 1st Qu.: NA
Median : 82.0 Median : 49.50 Median : 136.0 Median : NA
Mean :106.7 Mean : 64.34 Mean : 325.1 Mean :NaN
3rd Qu.:148.0 3rd Qu.: 77.75 3rd Qu.: 468.8 3rd Qu.: NA
Max. :286.0 Max. :329.00 Max. :1676.0 Max. : NA
NA's :1 NA's :38
MAULE ÑUBLE BBIO ARAUC RIOS
Min. : 5.0 Min. : NA Min. : 109.0 Min. : 1.0 Min. :203.0
1st Qu.: 99.0 1st Qu.: NA 1st Qu.: 253.5 1st Qu.: 1.0 1st Qu.:270.8
Median :127.0 Median : NA Median : 603.5 Median : 1.5 Median :319.5
Mean :165.3 Mean :NaN Mean : 848.0 Mean : 3.5 Mean :355.5
3rd Qu.:194.0 3rd Qu.: NA 3rd Qu.:1363.5 3rd Qu.: 4.0 3rd Qu.:473.8
Max. :526.0 Max. : NA Max. :2620.0 Max. :10.0 Max. :613.0
NA's :2 NA's :38 NA's :34 NA's :22
LAGOS AYSEN MAG ...17
Min. : 63.0 Min. : 1.00 Min. : 10.00 Min. : 803
1st Qu.: 661.2 1st Qu.: 8.00 1st Qu.: 28.25 1st Qu.:1753
Median : 842.0 Median :22.00 Median : 51.00 Median :3402
Mean :1267.2 Mean :20.36 Mean :112.82 Mean :3228
3rd Qu.:1799.5 3rd Qu.:29.50 3rd Qu.:157.75 3rd Qu.:4472
Max. :3515.0 Max. :41.00 Max. :411.00 Max. :6108
NA's :27 NA's :16
landing2 <- landing %>%
pivot_longer(cols=c("AyP" , "TPCA" ,"ANTOF", "ATCMA", "COQ",
"VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
"RIOS" , "LAGOS", "AYSEN", "MAG"),
names_to = "REGION",
values_to = "CAPTURA") %>%
rename(AÑO="...1") %>%
dplyr::select(-2)Genero un gráfico de barras
landing2$REGION <- factor(landing2$REGION ,
levels = c("AyP" , "TPCA" ,"ANTOF", "ATCMA", "COQ",
"VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
"RIOS" , "LAGOS", "AYSEN", "MAG"))
desem <- ggplot(landing2 %>%
drop_na(REGION),aes(AÑO, CAPTURA, fill=REGION)) +
geom_bar(stat="identity")+
scale_fill_viridis_d(option="G")+
theme_minimal()+
scale_x_continuous(breaks = seq(from = 1985, to = 2022, by = 4))+
scale_y_continuous(breaks = seq(from = 0, to = 4000, by = 1000))+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank(),
legend.position = "none")+
facet_wrap(~REGION, ncol=3)+
labs(y="Captura (t)",
x="")
desem
Los principales desembarques estan asociados a las regiones de Antifagasta, Valparaiso y BioBIo pero solo durante los primeros años (1985- mediados del 2000). Luego de esto, todas las regiones vieron disminuidos los registros de extracción del recurso.
BITÁCORA
Esta base tiene como principal objetivo obtener un indicador de esfuerzo para calcular un indice de abundancia relativo como la CPUE.
Primero identifico la estructura de la base y filtro el rrecurso bacalao
# filtro bitacoras dejando solo bacalao
bitb <- bit %>%
filter(COD_ESPECIE==37)%>%
mutate(dfp = as.double(dfp),
PESO_corr = as.double(PESO_corr),
prof_med=as.double(prof_med))
dim(bitb)
[1] 11074 77se eliminan 1208-11074 = 3134 registros
glimpse(bitb)ahora las estadísticas descriptivas de las variables de intéres. En este caso dfp, PESO_corr, prof_med
Miramos los oultiers de los datos
# OUTLIERS
#Order data
bitb <- bitb %>%
mutate(order = seq(1:nrow(bitb)))
#Select continuous variables to plot
p1 <- multi_dotplot(bitb, order, dfp)
p2 <- multi_dotplot(bitb, order, PESO_corr)
p3 <- multi_dotplot(bitb, order, prof_med)
p4 <- multi_dotplot(bitb, order, N_TRIPULANTES)
p5 <- multi_dotplot(bitb, order, HORAS_REPOSO)
p6 <- multi_dotplot(bitb, order, NUMERO_DE_ANZUELOS)
#Plot as a grid
ggarrange(p1, p2, p3, p4,p5, p6, nrow = 3)
Remuevo outliers
Los datos de profundidades andan bien, Filtro los datos de la variable dfp entre 0 y 30 días. y repito la inspección
bitb2 <- bitb %>%
filter(dfp>0,
dfp<80,
PESO_corr<50000,
N_TRIPULANTES<30,
NUMERO_DE_ANZUELOS<40000)Miro nuenvamente los datos filtrados
# OUTLIERS
#Order data
bitb2 <- bitb2 %>%
mutate(order = seq(1:nrow(bitb2)))
#Select continuous variables to plot
p1a <- multi_dotplot(bitb2, order, dfp)
p2a <- multi_dotplot(bitb2, order, PESO_corr)
p3a <- multi_dotplot(bitb2, order, prof_med)
p4a <- multi_dotplot(bitb2, order, N_TRIPULANTES)
p5a <- multi_dotplot(bitb2, order, HORAS_REPOSO)
p6a <- multi_dotplot(bitb2, order, NUMERO_DE_ANZUELOS)
#Plot as a grid
ggarrange(p1a, p2a, p3a, p4a, p5a, p6a, nrow = 3)
Identifico las dimensionesd e la nueva base y la comparo con la anterior
dim(bitb2)
[1] 4745 78
dim(bitb)
[1] 11074 78
table(bitb2$año)
2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
2 3 59 396 441 364 754 750 468 497 561 450Los registros se minimiza en ciertos años lo cual genera problemas en la fiabilidad de la estimación.
ahora calculamos el rendimiento nominal
CPUE Nominal
Ahora estimo la CPUE, la cual es el rendimoento entre el esfuerzo y la captura. La medida de esfuerzo será dfp. También calculo un rendimiento con la variable de esfuerzo HORAS_REPOSO pero solo como antecedente dado que no tengo muchos datos.
bitb2 <- bitb2 %>%
mutate(CPUE = PESO_corr/dfp) %>%
mutate(CPUE2 =PESO_corr/HORAS_REPOSO)
bitb2$prof_med_cat <- cut(bitb2$prof_med,
breaks = c(500, 1000, 1500,
2000, 3000))# Frequency polygon plot for catch
his <- ggplot(bitb2, aes(dfp)) +
geom_freqpoly(bins = 5) +
labs(x = "dfp", y = "Frequency") +
theme_bw() +
theme(panel.border = element_rect(colour = "black",
fill=NA, size = 1))
# Patterns in the variance? (any lack of homogeneity)
point <- ggplot(bitb2 %>%
filter(PESO_corr>1000), aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 5, alpha = 0.6) +
stat_smooth(method = "lm")+
theme(panel.background = element_blank()) +
theme(panel.border = element_rect(fill = NA, size = 1)) +
theme(strip.background = element_rect(fill = "white",
color = "white", size = 1)) +
theme(text = element_text(size=13)) +
xlab("dfp") + ylab("Peso Corregido")
#Plot as a grid
ggarrange(his, point, nrow = 1)
Tendencia del esfuerzo dfp a través de los años y por región
effor <- ggplot(bitb2 %>%
group_by(año, REGION_PUERTO_RECALADA) %>%
summarise(DFPM=mean(dfp)),
aes(año, DFPM))+
geom_point(shape = 16, size = 3, alpha = 0.7)+
scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 1))+
geom_smooth(method = 'lm',
colour = 'blue',
size = 1.5)+
facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
theme_few()+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank())+
labs(y="Effort (dfp)",
x="")
effor
Tendencia del esfuerzo NUMERO_DE_ANZUELOS a través de los años y por región
efforanz <- ggplot(bitb2 %>%
group_by(año, REGION_PUERTO_RECALADA) %>%
summarise(NANZ=mean(NUMERO_DE_ANZUELOS)),
aes(año, NANZ))+
geom_point(shape = 16, size = 3, alpha = 0.7)+
scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 1))+
geom_smooth(method = 'lm',
colour = 'green',
size = 1.5)+
facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
theme_few()+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank())+
labs(y="Effort (dfp)",
x="")
efforanz
Calculo los ceros de dfp
#CALCULATE NUMBER OF ZEROS
# What is the percentage of zeros i the response variable
round(sum(bitb2$dfp == 0) * 100 / nrow(bitb2),0)
[1] 0
#0Interacciones
# Interactions
# Year x season
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
theme_bw() +
xlab("Points") + ylab("Catch") +
facet_grid(año~trim)
# No
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
theme_bw() +
xlab("dfp") + ylab("Catch") +
facet_grid(año~N_TRIPULANTES)
# Perhaps
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
theme_bw() +
xlab("dfp") + ylab("Catch") +
facet_grid(trim~N_TRIPULANTES)
# CPUE slope varies between habitats - interactionProfundidad por region
profu <- ggplot(bitb2, aes(REGION_PUERTO_RECALADA, desc(prof_med),
group=REGION_PUERTO_RECALADA))+
geom_boxplot(fill=2, alpha=.5)+
geom_jitter(size=0.4, alpha=0.2,
width = .25)+
#facet_wrap(~año)+
theme_bw()+
scale_x_continuous(breaks = seq(from = 1, to = 15, by = 1))+
labs(x="Región",
y="Profundidad (mts)")
profu
ggplot(bitb2, aes(x = nombre_puerto_recalada , y = dfp, fill = nombre_puerto_recalada )) +
geom_violindot(fill_dots = "black") +
theme_modern() +
scale_fill_material_d()
ggplot(bitb2, aes(x = as.factor(REGION_PUERTO_RECALADA ),
y = bitb2$NUMERO_DE_ANZUELOS,
fill = as.factor(REGION_PUERTO_RECALADA))) +
geom_violin() +
theme_modern() +
scale_fill_material_d(palette="ice",
name="REGION")
ggplot(bitb2 %>%
drop_na(prof_med_cat),
aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm',
colour = 'red',
se = FALSE) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(prof_med_cat~año)
ploteo los datos de CPUE totales.
cpuenom <- ggplot(bitb2 %>%
group_by(zona, año) %>%
summarise(CPUEM=mean(CPUE)) %>%
filter(CPUEM<600,
año>1996,
zona!=4),
aes(año, CPUEM, color=zona, group=zona))+
geom_point()+
geom_smooth()+
scale_color_viridis_b(option="B",
labels = seq(1, 3, by = 1))+
scale_x_continuous(breaks = seq(from = 2010, to = 2022, by = 1))+
theme_few()+
labs(x="",
y="CPUE")
cpuenom
guardo los datos bitacora CPUE nominal 1996-2022 del objeto cpuenom
write_csv(cpuenom, "CPUE_MON.csv")CPUE Estandarizada
Identificar los proncipales factores para modelar la variable.
Luego de una reunión con Seguimiento (Patricio Galvez) se indicaron algunos aspectos que deben ser considerados en la estandariación.
El poder de pesca (tipo de embarcación) es mayor en la zona centro sur de Chile. Identificar equilibrio de la base
bitb2respecto a este dato.Usar zonas como factor en desmedro de región.
Identifico una estadistica descriptiva general de la base de estandarización
Correlaciones
Análisis de utilidad para identificar a traves de un metodo de corrrelación de pearson, la correlación en tre las variables que serán utilizadas en la estandariacion de la CPUE.
Primero identifico las variables a correlacionar
names(bitb2)
[1] "ID" "embarque"
[3] "COD_BARCO" "FECHA_HORA_RECALADA"
[5] "FECHA_HORA_ZARPE" "COD_PESQUERIA"
[7] "COD_PESQUERIA_ant" "N_TRIPULANTES"
[9] "N_CALADAS" "PUERTO_ZARPE"
[11] "PUERTO_RECALADA" "nombre_puerto_recalada"
[13] "NRO_FORMULARIO" "REGION_PUERTO_RECALADA"
[15] "ESPECIE_OBJETIVO" "ECOSONDA"
[17] "GPS" "POTENCIA_MOTOR_EQ"
[19] "VIRADOR" "GASTOS_VIAJE"
[21] "NUMERO_LANCE_EX" "FECHA_LANCE"
[23] "año" "mes"
[25] "trim" "dfp"
[27] "ID_CUADRICULA" "ID_PROCEDENCIA"
[29] "PESO_TOTAL_CAPTURA" "LATITUD"
[31] "LONGITUD" "lat_ctgr"
[33] "lon_ctgr" "zona"
[35] "ESPECIE_OBJETIVO_LANCE" "CLASE_LANCE"
[37] "COD_ESPECIE" "PESO"
[39] "PESO_corr" "CAPTURA_RETENIDA"
[41] "VOLUMEN_DESCARTE" "LUGAR_DESCARTE"
[43] "TIPO_DESCARTE" "PRECIO_UNITARIO"
[45] "DESTINO_CAPTURA" "NRO_ANZUELOS"
[47] "REGION_PROCEDENCIA" "FECHA_HORA_FIN_CALADO"
[49] "FECHA_HORA_INI_CALADO" "LATITUD_INICIO_CALADO"
[51] "LONGITUD_INICIO_CALADO" "LATITUD_FIN_CALADO"
[53] "LONGITUD_FIN_CALADO" "FECHA_HORA_INI_VIRADO"
[55] "FECHA_HORA_FIN_VIRADO" "LATITUD_INICIAL_VIRADO"
[57] "LONGITUD_INICIAL_VIRADO" "LATITUD_FIN_VIRADO"
[59] "LONGITUD_FIN_VIRADO" "PROFUNDIDAD_MINIMA_LM"
[61] "PROFUNDIDAD_MAXIMA_LM" "PROFUNDIDAD_FONDO_VIRADO_INI"
[63] "PROFUNDIDAD_FONDO_VIRADO_FIN" "prof_med"
[65] "SEPARACION_ANZUELO" "MARCA_ANZUELO"
[67] "NUMERO_DE_ANZUELOS" "PERDIDA_ANZUELOS"
[69] "TIPO_CARNADA" "TAMANIO_ANZUELO"
[71] "INTENSIDAD_VIENTO" "LONGITUD_LINEA_MADRE"
[73] "HORAS_REPOSO" "TIPO_ESPINEL"
[75] "N_PANOS" "ANZ_POR_PANO"
[77] "filtro" "order"
[79] "CPUE" "CPUE2"
[81] "prof_med_cat"
bitb3 <- bitb2 %>%
dplyr::select(8, 18, 23,24, 25, 26, 32, 33, 34, 39, 69, 73, 75, 78, 79, 80, 81)results <- correlation(bitb3)
results %>%
summary(redundant=TRUE) %>%
plot(results, show_data = "points")+
theme_bw()
skimr::skim(bitb3)| Name | bitb3 |
| Number of rows | 4745 |
| Number of columns | 17 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| factor | 1 |
| numeric | 15 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| TIPO_CARNADA | 122 | 0.97 | 1 | 5 | 0 | 8 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| prof_med_cat | 939 | 0.8 | FALSE | 4 | (1e: 1867, (50: 961, (1.: 853, (2e: 125 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| N_TRIPULANTES | 0 | 1.00 | 7.23 | 1.62 | 2.00 | 6.00 | 7.00 | 9.00 | 10.00 | ▁▂▆▇▆ |
| POTENCIA_MOTOR_EQ | 167 | 0.96 | 405.74 | 399.03 | 0.00 | 320.00 | 390.00 | 420.00 | 3903.00 | ▇▁▁▁▁ |
| año | 0 | 1.00 | 2018.05 | 2.47 | 2011.00 | 2016.00 | 2018.00 | 2020.00 | 2022.00 | ▁▅▆▆▇ |
| mes | 0 | 1.00 | 7.29 | 3.79 | 1.00 | 3.00 | 9.00 | 11.00 | 12.00 | ▅▂▁▂▇ |
| trim | 0 | 1.00 | 2.71 | 1.32 | 1.00 | 1.00 | 3.00 | 4.00 | 4.00 | ▆▂▁▂▇ |
| dfp | 0 | 1.00 | 16.73 | 6.62 | 0.67 | 13.02 | 16.79 | 19.79 | 50.67 | ▂▇▂▁▁ |
| lat_ctgr | 51 | 0.99 | -35.49 | 6.89 | -55.17 | -40.99 | -35.25 | -33.17 | -18.33 | ▁▅▇▂▂ |
| lon_ctgr | 51 | 0.99 | -73.24 | 1.62 | -81.89 | -74.69 | -73.08 | -71.97 | -51.94 | ▁▇▁▁▁ |
| zona | 0 | 1.00 | 2.11 | 0.65 | 1.00 | 2.00 | 2.00 | 2.00 | 4.00 | ▂▇▁▃▁ |
| PESO_corr | 0 | 1.00 | 755.69 | 1304.56 | 0.00 | 41.00 | 110.00 | 1079.66 | 14092.00 | ▇▁▁▁▁ |
| HORAS_REPOSO | 2988 | 0.37 | 1629.72 | 832.08 | 12.00 | 1200.00 | 1200.00 | 2400.00 | 7200.00 | ▇▅▁▁▁ |
| N_PANOS | 1833 | 0.61 | 226.33 | 280.51 | 1.00 | 1.00 | 2.00 | 420.00 | 900.00 | ▇▁▃▁▁ |
| order | 0 | 1.00 | 2373.00 | 1369.91 | 1.00 | 1187.00 | 2373.00 | 3559.00 | 4745.00 | ▇▇▇▇▇ |
| CPUE | 0 | 1.00 | 54.95 | 92.65 | 0.00 | 2.44 | 6.49 | 84.01 | 1067.68 | ▇▁▁▁▁ |
| CPUE2 | 2988 | 0.37 | 1.55 | 5.23 | 0.00 | 0.51 | 0.99 | 1.79 | 158.44 | ▇▁▁▁▁ |
Datos 0
Are there missing values? y sirve para identificar los factores que se utilzaran en la estandarización.
dim(bitb3)
[1] 4745 17
colSums(is.na(bitb3))
N_TRIPULANTES POTENCIA_MOTOR_EQ año mes
0 167 0 0
trim dfp lat_ctgr lon_ctgr
0 0 51 51
zona PESO_corr TIPO_CARNADA HORAS_REPOSO
0 0 122 2988
N_PANOS order CPUE CPUE2
1833 0 0 2988
prof_med_cat
939Es probable que HORAS_REPOSOde la red no podamos usarlo como factor ni como variable de esfuerzo por lo exiguo del dato.
Variabilidad
asd1<-ggplot(data = bitb3 %>%
filter(zona!=4),
aes(y=dfp, x=zona,
group=zona)) +
geom_boxplot()+
labs(title = "Variabilidad en DFP",
y = "DFP", x = "ZONA") +
facet_wrap(~ año, scales = "free")+
theme_few()
asd1
asd2<-ggplot(data =bitb3%>%
drop_na(prof_med_cat) %>%
filter(zona!=4),
aes(x=prof_med_cat,
y=lat_ctgr)) +
geom_boxplot()+
labs(title = "Variabilidad en distribución latitudinal",
y = "Latitud", x = "Prof") +
theme_few()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
facet_wrap(~ año)
asd2
Normality and Homogenity factors
# Frequency polygon plot for catch
ggplot(bitb3 %>%
filter(zona!=4),
aes(CPUE)) +
geom_freqpoly(bins = 8) +
labs(x = "Cpue", y = "Frequency") +
facet_wrap(.~zona)+
theme_few()+
theme(panel.border =
element_rect(colour = "black",
fill=NA, size = 1))
Shapiro-Wilk test for deviation from normality
shapiro.test(bitb3$CPUE)
Shapiro-Wilk normality test
data: bitb3$CPUE
W = 0.63039, p-value < 2.2e-16Departure from normality…but not critical
Colinealidad entre factores
# COLLINEARITY
Coll <- c(
"año", "dfp" ,"lat_ctgr", "zona" ,
"N_PANOS" , "CPUE", "prof_med_cat")
# Obtain summary using the ggpairs command from the GGally library
ggpairs(bitb3[,Coll], ggplot2::aes(alpha = 0.9, colour=prof_med_cat),
lower = list(combo = "count"))+
scale_fill_manual(values=c("red", "blue", "green", "black")) +
scale_colour_manual(values=c("red", "blue", "green", "black")) +
theme_few()
# Nothing serious#Calculate Variance Inflation Factor (VIF)
vifmodel <- round(vif(lm(PESO_corr ~ zona + dfp + año + mes,
data = bitb3)),2)
vifmodel
zona dfp año mes
1.34 1.29 1.06 1.01barplot(vifmodel, main = "VIF Values", horiz = TRUE, col = "steelblue")
La multicolinealidad en el análisis de regresión se produce cuando dos o más variables predictoras están altamente correlacionadas entre sí, de modo que no proporcionan información única o independiente en el modelo de regresión.
Si el grado de correlación es lo suficientemente alto entre las variables, puede causar problemas al ajustar e interpretar el modelo de regresión.
La forma más común de detectar la multicolinealidad es mediante el factor de inflación de varianza (VIF), que mide la correlación y la fuerza de la correlación entre las variables predictoras en un modelo de regresión.
El valor de VIF comienza en 1 y no tiene límite superior. Una regla general para interpretar VIF es la siguiente:
Un valor de 1 indica que no hay correlación entre una variable predictora dada y cualquier otra variable predictora en el modelo. Un valor entre 1 y 5 indica una correlación moderada entre una variable predictora dada y otras variables predictoras en el modelo, pero esto a menudo no es lo suficientemente grave como para requerir atención. Un valor superior a 5 indica una correlación potencialmente grave entre una variable predictora dada y otras variables predictoras en el modelo. En este caso, es probable que las estimaciones de coeficientes y los valores p en el resultado de la regresión no sean fiables.
Interaccion
# Patterns in the variance? (Evidence for lack of homogeneity)
ggplot(bitb3,
aes(x = dfp,
y = (CPUE))) +
geom_jitter(shape = 16,
size = 2.5,
alpha = 0.3,
height = 0.25,
width = 0.5) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
facet_grid(bitb3$N_TRIPULANTES~bitb3$zona)+
theme(panel.background = element_blank()) +
theme(panel.border = element_rect(fill = NA,
size = 1)) +
theme(strip.background =
element_rect(fill = "white",
color = "white",
size = 1)) +
theme(text = element_text(size=13)) +
theme_few()+
xlab("Effort") +
ylab("CPUE") +
ggtitle("variabless por Numero de Tripulantes y por zona")
Year x mes
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(año~trim)
# No# Year x habitat
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(zona~prof_med_cat)
# Perhaps# Zona
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(año~trim)
facet_grid(~zona)
<ggproto object: Class FacetGrid, Facet, gg>
compute_layout: function
draw_back: function
draw_front: function
draw_labels: function
draw_panels: function
finish_data: function
init_scales: function
map_data: function
params: list
setup_data: function
setup_params: function
shrink: TRUE
train_scales: function
vars: function
super: <ggproto object: Class FacetGrid, Facet, gg>
# CPUE slope varies between habitats - interaction# Tipo carnada zona
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few()+
xlab("Points") + ylab("Catch") +
facet_grid(TIPO_CARNADA~zona)
# PerhapsGenero categorias y Factorizo
bitb3 <- bitb3 %>%
mutate(POTENCIA_CAT = cut(POTENCIA_MOTOR_EQ,
breaks = c(0, 200, 300, 400, 1000),
labels = c("Bajo", "Medio-Bajo", "Medio-Alto", "Alto")),
LAT_CAT = cut(lat_ctgr,
breaks = c(-20, -35, -40, -47),
labels = c("Norte", "Centro", "Sur")),
N_TRIPULANTES=factor(N_TRIPULANTES),
año=factor(año),
zona=factor(zona),
trim=factor(trim)) %>%
drop_na(zona,
LAT_CAT,
prof_med_cat,
POTENCIA_CAT) %>%
filter(CPUE>0) %>%
mutate(logCPUE=log(CPUE))
glimpse(bitb3)
Rows: 3,347
Columns: 20
$ N_TRIPULANTES <fct> 4, 8, 8, 8, 6, 7, 6, 7, 8, 8, 6, 8, 8, 9, 9, 9, 9, 9…
$ POTENCIA_MOTOR_EQ <dbl> 250, 380, 320, 450, 380, 420, 370, 240, 320, 360, 45…
$ año <fct> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
$ mes <dbl> 7, 7, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11,…
$ trim <fct> 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
$ dfp <dbl> 13.833333, 11.741667, 15.086806, 13.729167, 18.31736…
$ lat_ctgr <dbl> -40.66667, -39.96667, -40.50000, -33.91667, -40.6666…
$ lon_ctgr <dbl> -74.50000, -74.33333, -74.50000, -72.58333, -74.6666…
$ zona <fct> 2, 2, 2, 2, 2, 3, 3, 3, 3, 2, 3, 3, 3, 2, 2, 2, 2, 2…
$ PESO_corr <dbl> 1095.0568, 1256.6812, 1422.8584, 975.6000, 2026.6464…
$ TIPO_CARNADA <chr> "2", "3", "3", NA, "1", "2", "2", "2", "3", "3", "3"…
$ HORAS_REPOSO <dbl> 1400, 1800, 1200, 4800, 1200, 1200, 2400, 1200, 1200…
$ N_PANOS <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ order <int> 8, 10, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 2…
$ CPUE <dbl> 79.160730, 107.027493, 94.311443, 71.060393, 110.640…
$ CPUE2 <dbl> 0.7821834, 0.6981562, 1.1857153, 0.2032500, 1.688872…
$ prof_med_cat <fct> "(1e+03,1.5e+03]", "(500,1e+03]", "(1e+03,1.5e+03]",…
$ POTENCIA_CAT <fct> Medio-Bajo, Medio-Alto, Medio-Alto, Alto, Medio-Alto…
$ LAT_CAT <fct> Norte, Centro, Norte, Sur, Norte, Norte, Norte, Nort…
$ logCPUE <dbl> 4.371480, 4.673086, 4.546603, 4.263530, 4.706288, 4.…The data exploration showed:
- One significant outlier in catch
- Deviation from normality in response variable
- Possible departure from homogeneity
- Few zeros in the response variable
- No collinearity
- No imbalance
- Possible season x habitat interaction
Chequeo distribuciones
bcpue <- ggplot(bitb3,
aes(CPUE))+
coord_flip()+
geom_boxplot()+
theme_few()
hcpue <- ggplot(bitb3,
aes(CPUE))+
geom_histogram(bins=15, fill=2)+
theme_few()
ggarrange(bcpue, hcpue, ncol = 2) 
blcpue <- ggplot(bitb3,
aes(logCPUE))+
coord_flip()+
geom_boxplot()+
theme_few()
hlcpue <- ggplot(bitb3,
aes(logCPUE))+
geom_histogram(bins=15, fill=2)+
theme_few()
ggarrange(blcpue, hlcpue, ncol = 2) 
#Chequeo visual de normalidad y dispersion
#-------------------------------------------------------------------------------
#CPUE y LOGCPUE
par(mfrow=c(2,2))
qqnorm(bitb3$CPUE)
qqline(bitb3$CPUE,col='red') #
boxplot(CPUE ~ año,data=bitb3,na.rm=T,#
main="Boxplot nominal CPUE by year", ylab="CPUE",col='light grey',boxwex=0.65)
qqnorm(bitb3$logCPUE)
qqline(bitb3$logCPUE,col='red') #
boxplot(logCPUE ~ año,data=bitb3,na.rm=T,#
main="Boxplot nominal LogCPUE by year", ylab="LogCPUE",col='light grey',boxwex=0.65)
Modelos GLM
#Predictores (Factores Principales)
gauss01 = glm(formula= logCPUE ~ año,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss02 = glm(formula= logCPUE ~ año+trim,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss03 = glm(formula= logCPUE ~ año+trim+prof_med_cat,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss04 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss05 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss06 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+TIPO_CARNADA,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss07 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)
gauss08 = glm(formula= logCPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona:prof_med_cat,
family = gaussian(link="identity"),
data = bitb3,na.action=na.exclude)rmodelo01 = c(AIC(gauss01),(gauss01$null.deviance-gauss01$deviance)/gauss01$null.deviance)
rmodelo02 = c(AIC(gauss02),(gauss02$null.deviance-gauss02$deviance)/gauss02$null.deviance)
rmodelo03 = c(AIC(gauss03),(gauss03$null.deviance-gauss03$deviance)/gauss03$null.deviance)
rmodelo04 = c(AIC(gauss04),(gauss04$null.deviance-gauss04$deviance)/gauss04$null.deviance)
rmodelo05 = c(AIC(gauss05),(gauss05$null.deviance-gauss05$deviance)/gauss05$null.deviance)
rmodelo06 = c(AIC(gauss06),(gauss06$null.deviance-gauss06$deviance)/gauss06$null.deviance)
rmodelo07 = c(AIC(gauss07),(gauss07$null.deviance-gauss07$deviance)/gauss07$null.deviance)
rmodelo08 = c(AIC(gauss08),(gauss08$null.deviance-gauss08$deviance)/gauss08$null.deviance)
resultados = rbind(rmodelo01,rmodelo02,rmodelo03,rmodelo04,rmodelo05,rmodelo06,
rmodelo07,rmodelo08);resultadoscheck_model(gauss07)
Table with significance level (***).
tab_model(gauss01,
gauss02,
gauss03,
gauss04,
gauss05,
gauss06,
gauss07,
gauss08,
p.style = "stars")model_performance(gauss07)
# Indices of model performance
AIC | AICc | BIC | R2 | RMSE | Sigma
---------------------------------------------------------
12846.652 | 12847.141 | 13017.894 | 0.335 | 1.635 | 1.642Table comparing performance model.
compare_performance(gauss01,
gauss02,
gauss03,
gauss04,
gauss05,
gauss06,
gauss07,
gauss08,
rank = TRUE,
verbose = FALSE)
# Comparison of Model Performance Indices
Name | Model | R2 | RMSE | Sigma | Performance-Score
-----------------------------------------------------------
gauss08 | glm | 0.341 | 1.627 | 1.635 | 100.00%
gauss07 | glm | 0.335 | 1.635 | 1.642 | 97.35%
gauss06 | glm | 0.334 | 1.643 | 1.650 | 95.42%
gauss05 | glm | 0.299 | 1.679 | 1.685 | 82.38%
gauss04 | glm | 0.264 | 1.720 | 1.726 | 68.13%
gauss03 | glm | 0.217 | 1.774 | 1.778 | 49.32%
gauss02 | glm | 0.102 | 1.900 | 1.904 | 4.17%
gauss01 | glm | 0.091 | 1.912 | 1.915 | 0.00%Regarding this result, we can select best model and check some assumtion like:
- Colinealidad
- Data points influyentes
- Homoscedasticidad
- Normalidad de los residuos
- Independencia
COMPOSICIONES DE TALLAS
La información contenida en las estructiras de tallas son consideradas una de las piezas mas importantes en la ciencia pesquera (Canales, Punt, & Mardones, 2021; Hordyk, Ono, Sainsbury, Loneragan, & Prince, 2014; Rudd & Thorson, 2018). Sin embargo, esta debe ser validada y confirmada para no violar supuestos referidos a la repreentatividad del dato y su relación con la dinámica poblacional.
Ahora nos disponemos a explorar esta fuente de información.
Primero identificamos la estructura de la base;
glimpse(long)
Rows: 15,831
Columns: 24
$ ID <chr> "COD_BARCO=730881 AND TO_CHAR(FECHA_HORA_RECALADA,…
$ COD_BARCO <dbl> 730881, 730881, 730881, 730881, 730881, 730881, 73…
$ FECHA_HORA_RECALADA <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
$ FECHA_HORA_ZARPE <dttm> 2021-05-10 11:00:00, 2021-05-10 11:00:00, 2021-05…
$ COD_PESQUERIA <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50…
$ PUERTO_RECALADA <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 3,…
$ NRO_FORMULARIO <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,…
$ REGION <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1,…
$ FECHA_LANCE <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
$ NUMERO_LANCE_EX <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ LATITUD <dbl> 243900, 243900, 243900, 243900, 243900, 243900, 24…
$ LONGITUD <dbl> 714400, 714400, 714400, 714400, 714400, 714400, 71…
$ NUMERO_EJEMPLARES <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PESO_TOTAL_MUESTRA <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ NUMERO_CAJA <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ COD_ESPECIE <dbl> 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37…
$ ORIGEN_MUESTRA <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3,…
$ FECHA_MUESTREO <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ N_TOTAL_INDIV <dbl> 85, 85, 85, 85, 85, 85, 85, 84, 84, 84, 84, 84, 23…
$ LONGITUD_DESCARTE <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LONGITUD_MUESTRA <dbl> 115, 125, 129, 130, 143, 152, 157, 100, 108, 110, …
$ SEXO <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ N_INDIVIDUOS <dbl> 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 1, 1,…
$ año <dbl> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 20…
colSums(is.na(long))
ID COD_BARCO FECHA_HORA_RECALADA FECHA_HORA_ZARPE
2997 0 0 0
COD_PESQUERIA PUERTO_RECALADA NRO_FORMULARIO REGION
0 0 0 0
FECHA_LANCE NUMERO_LANCE_EX LATITUD LONGITUD
0 0 4425 4425
NUMERO_EJEMPLARES PESO_TOTAL_MUESTRA NUMERO_CAJA COD_ESPECIE
15733 15831 0 0
ORIGEN_MUESTRA FECHA_MUESTREO N_TOTAL_INDIV LONGITUD_DESCARTE
0 15831 0 15831
LONGITUD_MUESTRA SEXO N_INDIVIDUOS año
0 0 2 0Primero selecciono las columnas de interes y luego genero la expansión de LONGITUD_MUESTRA a N_INDIVIDUOS. Expand frecuency data related length, in this case N_INDIVIDUOS column have frecuency that we need expand to whole data frame.
Selecciono variables de interes
long1 <- long %>%
dplyr::select(6,8,11,12,21,22,23,24)
dim(long1)
[1] 15831 8long2 <- long %>%
drop_na(N_INDIVIDUOS) %>%
type.convert(as.is = TRUE) %>%
uncount(N_INDIVIDUOS)
dim(long2)
[1] 34299 23
legend.labels <- c('Indet.', 'Macho' , 'Hembra')
nbco <- ggplot(long2, aes(x=LONGITUD_MUESTRA, y = as.factor(año),
fill=as.factor(SEXO)))+
geom_density_ridges(stat = "density_ridges", bins = 20,
scale = 2, draw_baseline = FALSE,
alpha=0.5)+
facet_wrap(.~REGION, ncol=5) +
geom_vline(xintercept = 110, color = "red")+
scale_fill_viridis_d(option = "C",
name="SEXO",
labels=legend.labels)+
scale_y_discrete(breaks = seq(from = 2004, to = 2022, by = 2))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
theme_minimal()+
xlab("Longitud (cm.)")+
ylab("")
#scale_x_discrete((limits = rev(levels(talla2021$ANO_ARR))))+
nbco
Prepara los vectores para sumar a los .dat del modelo si los necesito.
long2$LONG_CAT <- as.numeric(as.character(cut(x = long2$LONGITUD_MUESTRA,
breaks = seq(10,220,2),
labels = seq(10,218,2),
right = FALSE)))
LONGT <- table(long2$año, long2$LONG_CAT)
LONGT
10 12 38 40 42 44 48 50 52 54 56 58 60 62 64 66 68 70
2004 2 1 0 0 0 0 0 0 0 0 0 3 14 17 3 13 14 23
2005 0 0 0 0 0 0 0 0 0 0 1 0 1 2 5 10 11 17
2006 0 0 0 0 0 0 0 0 0 2 2 3 2 4 6 6 8 12
2007 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2008 0 0 0 0 0 0 0 0 1 0 0 0 2 3 2 3 10 3
2009 0 0 0 0 0 0 0 0 0 0 1 0 5 2 4 8 25 23
2010 0 0 0 0 0 0 0 0 0 1 1 2 1 5 7 11 11 16
2011 0 0 0 1 2 1 1 1 0 1 2 2 4 9 12 8 17 21
2012 0 0 0 0 0 0 0 0 0 0 0 0 1 5 2 20 20 44
2013 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 12 10
2014 0 0 1 0 0 0 1 0 0 0 1 2 27 21 41 85 112 120
2015 0 0 0 0 0 0 0 0 0 0 0 1 6 16 20 41 65 87
2016 0 0 0 0 0 0 0 0 0 0 1 4 10 17 47 56 90 140
2017 0 0 1 0 0 0 0 0 0 3 9 19 60 59 122 117 111 178
2018 0 0 0 0 0 0 0 0 0 6 8 25 27 35 46 61 47 68
2019 0 0 0 0 0 0 0 0 0 1 2 4 18 42 81 118 139 202
2020 0 0 0 0 0 0 0 0 0 0 0 0 4 5 11 17 11 20
2021 0 0 0 0 0 0 0 0 1 0 0 0 3 4 10 10 18 41
2022 0 0 0 0 0 0 0 0 2 1 5 13 24 32 48 83 88 139
72 74 76 78 80 82 84 86 88 90 92 94 96 98 100 102 104 106
2004 22 22 40 28 33 50 31 62 54 72 47 37 41 43 52 21 17 22
2005 19 19 26 21 21 31 31 22 31 25 22 18 15 22 10 17 9 8
2006 11 13 17 20 26 13 12 11 9 9 7 3 12 6 6 9 8 7
2007 0 1 0 0 1 1 1 2 0 2 0 0 2 0 0 0 0 0
2008 17 9 17 22 12 15 13 17 16 10 12 10 9 10 6 7 6 10
2009 39 52 83 82 84 113 105 94 97 82 98 71 54 50 38 28 43 15
2010 14 22 31 36 21 34 36 41 31 57 36 28 28 20 30 20 11 15
2011 11 20 25 21 11 18 13 12 18 9 12 10 13 13 8 9 5 8
2012 63 67 60 39 61 32 29 24 18 27 34 14 12 6 13 9 14 5
2013 29 29 37 61 99 82 95 112 62 115 55 70 68 41 44 42 34 43
2014 186 188 216 230 240 252 243 241 224 190 166 124 103 101 68 70 51 36
2015 116 165 186 186 235 248 253 268 195 240 237 193 179 142 139 86 81 47
2016 147 192 176 223 278 280 270 270 228 239 234 228 166 154 154 140 99 88
2017 159 171 146 152 158 160 196 178 147 190 162 184 180 160 217 131 139 94
2018 46 50 33 28 35 25 23 43 33 43 32 45 44 41 37 36 30 41
2019 217 190 165 124 102 80 54 31 28 39 24 33 26 8 43 23 26 27
2020 29 30 28 19 21 23 13 1 6 8 3 4 1 7 4 4 2 7
2021 54 77 104 115 179 143 159 128 107 114 75 61 48 30 36 25 19 14
2022 173 203 224 270 323 372 383 376 303 335 285 252 192 119 123 119 76 56
108 110 112 114 116 118 120 122 124 126 128 130 132 134 136 138 140 142
2004 15 16 11 9 14 11 21 7 6 2 4 8 9 5 4 2 6 3
2005 12 10 8 7 5 3 3 3 2 5 1 2 1 0 2 1 2 0
2006 6 2 5 5 2 2 4 4 3 2 1 0 1 3 1 0 0 2
2007 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0
2008 7 6 7 11 6 6 2 4 2 4 5 4 2 2 3 1 3 3
2009 23 25 18 19 15 12 13 9 12 8 9 10 4 8 11 10 15 11
2010 16 13 6 8 5 7 9 4 4 5 0 1 1 1 2 4 1 2
2011 7 15 15 4 1 7 2 4 9 4 4 2 1 1 2 3 2 1
2012 8 10 3 3 2 3 2 0 1 1 0 0 1 0 1 0 2 0
2013 31 35 21 21 15 12 21 8 7 7 6 9 5 4 3 4 5 0
2014 46 35 38 38 41 35 41 29 26 13 9 17 12 8 9 7 9 11
2015 52 52 39 34 44 41 33 19 27 29 16 18 18 8 15 12 15 5
2016 45 58 42 30 30 29 28 22 20 23 11 18 13 9 5 6 7 5
2017 77 139 60 78 45 37 55 19 24 20 11 26 18 22 17 5 12 10
2018 30 25 26 22 14 10 21 14 4 9 6 10 2 9 9 5 5 7
2019 24 26 20 22 7 17 22 15 14 11 7 5 1 7 8 1 6 3
2020 3 3 2 5 4 4 4 0 2 0 4 2 1 0 1 1 1 0
2021 13 26 34 20 29 21 32 21 27 27 20 17 21 12 20 15 17 13
2022 33 52 22 21 35 29 14 22 26 21 37 24 21 27 25 14 19 20
144 146 148 150 152 154 156 158 160 162 164 166 168 170 172 174 176 178
2004 0 3 2 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0
2005 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2006 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0
2007 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2008 1 3 5 1 3 1 3 3 0 3 0 0 0 0 0 0 0 0
2009 4 7 9 13 10 2 8 10 5 5 2 5 5 0 1 1 2 0
2010 5 0 2 3 0 2 1 2 4 0 0 0 1 1 1 0 0 0
2011 1 0 2 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0
2012 1 0 0 1 0 1 0 0 0 0 0 1 0 0 2 0 0 1
2013 5 4 2 0 0 2 1 0 0 0 2 1 0 0 0 1 0 0
2014 5 5 6 1 1 2 1 1 3 2 2 3 1 1 0 0 1 0
2015 10 6 10 9 5 5 3 4 2 2 0 3 3 3 2 1 1 0
2016 7 7 4 4 4 5 0 1 0 4 3 5 1 2 0 0 3 0
2017 8 9 7 9 8 4 0 5 6 6 4 5 0 3 2 1 2 0
2018 5 1 4 6 2 1 1 5 3 1 1 0 1 2 1 3 0 0
2019 3 4 2 4 1 1 3 0 2 0 2 0 0 0 0 0 0 0
2020 2 0 0 0 1 0 0 0 2 0 0 0 0 2 0 0 0 0
2021 11 7 5 13 7 2 5 5 8 6 2 2 4 5 0 0 0 1
2022 19 18 16 22 17 5 10 7 13 9 6 5 7 5 3 4 3 2
180 182 184 186 188 190 192 200 202
2004 1 0 0 0 0 0 0 0 0
2005 0 0 0 0 0 0 0 0 0
2006 0 0 0 0 0 0 0 0 0
2007 0 0 0 0 0 0 0 0 0
2008 0 0 0 0 0 0 0 0 0
2009 0 0 1 0 0 0 0 0 0
2010 0 1 0 0 0 0 0 0 0
2011 0 0 0 0 0 0 0 0 0
2012 0 0 0 0 0 0 0 0 0
2013 0 0 0 0 0 0 0 0 0
2014 0 0 0 0 0 0 1 0 0
2015 0 0 0 0 0 0 0 0 0
2016 0 0 1 0 1 0 0 0 0
2017 0 1 0 0 0 0 0 1 0
2018 0 0 0 0 0 0 0 0 1
2019 0 0 0 0 0 0 0 0 0
2020 0 0 0 0 0 0 0 0 0
2021 0 0 0 0 0 0 0 0 0
2022 0 0 3 1 0 1 0 0 0write.csv()BITACORA FIPA 96- 14
bitfipa<- bit8696 %>%
mutate(CPUE = CAPTURA/DFP) %>%
drop_na(CPUE)
cpuefipa <- ggplot(bitfipa %>%
filter(ANO<1997) %>%
group_by(ANO) %>%
summarise(CPUEF=mean(CPUE)),
aes(ANO, CPUEF))+
geom_point(shape = 16, size = 3, alpha = 0.7)+
scale_x_continuous(breaks = seq(from = 1986, to = 2022, by = 1))+
geom_smooth(method = 'lm',
colour = 'blue',
size = 1.5)+
# stat_regline_equation(label.x=1993, label.y=500)+
# stat_cor(aes(label=..rr.label..), label.x=1993, label.y=450)+
#facet_wrap(~AREA)+
theme_bw()+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank())+
labs(y="CPUE (kg/dfp)",
x="")
cpuefipa
Guardo datos FIPA 96
cpuenomfipa <- bitfipa %>%
filter(ANO<1997) %>%
group_by(ANO) %>%
summarise(CPUEF=mean(CPUE))write_csv(cpuenomfipa, "CPUE_FIPA.csv")MAPAS
Lo primero es transformar los datos de lon_ctgr y lat_ctgr a formarto geom. A su vez, saco las NA de las coord
transformar los datos en un sf object
codmap <- st_as_sf(bitb2 %>%
drop_na(lon_ctgr) %>%
drop_na(lat_ctgr),
coords = c("lon_ctgr", "lat_ctgr"), crs = 4326)Ahora genero los mapas de Chile con un raster.
chile <- raster::getData("GADM", country = "CHL", level = 0)
chile1<-fortify(chile)
chilemap <- ggplot()+
geom_polygon(data=chile, aes(x=long, y=lat, group=group),
fill="lightblue",color="grey20", size=0.15)+
coord_sf(crs = st_crs(4326),
xlim = c(-80, -65),
ylim = c(-58, -15)) +
theme_bw()
chilemap
# Aca veo los nombres
chile@data$NAME_1
NULLLuego genero los bordes sobre los cuales haré la grilla.
e <- extent(-80,-65,-58,-15)
rc <- crop(chile, e)
# la proyección adecuada en lat long
tcrs <- CRS("+init=epsg:4326")
transformed_points <- spTransform(rc, tcrs)
# para dejarlo en formato geom_sf
chile2 <- st_as_sf(transformed_points) Se estructura la grilla en el objeto raster de chile1
grid <- chile2 %>%
st_make_grid(cellsize = c(1,1)) %>% # para que quede cuadrada
st_cast("MULTIPOLYGON") %>%
st_sf() %>% # objeto en spatial feature
mutate(cellid = row_number())Pongo los datos codmap en la grilla. Aca solo elegí dfpy CPUE pero se ùeden resumir otros
joindat <- grid %>%
st_join(codmap) %>%
group_by(cellid, año) %>%
summarise(CPUEM = mean(CPUE),
DFPM =mean(dfp)) # por ejemplo, plotear la captura "CAPTURA_1"Mapa Esfuerzo
here we can viz diferent variables
## Plot final
mas1 <- ggplot() +
geom_sf(data=joindat %>%
filter(!is.na(DFPM)), aes(fill = DFPM),
color=NA) +
scale_fill_viridis_b(option="E",
direction=-1, name="DFP")+
geom_sf(data = grid, fill=NA, color=NA) +
geom_sf(data = chile2, color="grey", fill="white") +
coord_sf() +
geom_hline(yintercept = -47, color = "red")+
scale_alpha(guide="none")+
facet_wrap(~año, ncol=6)+
scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
guides(colour = guide_legend()) +
theme_bw()+
theme(panel.background = element_rect(fill = 'white'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
mas1
Mapa CPUE
## Plot final
mas2 <- ggplot() +
geom_sf(data=joindat %>%
filter(!is.na(CPUEM)), aes(fill = CPUEM),
color=NA) +
scale_fill_viridis_b(option="G",
direction=-1, name="CPUE")+
geom_sf(data = grid, fill=NA, color=NA) +
geom_sf(data = chile2, color="grey", fill="white") +
coord_sf() +
geom_hline(yintercept = -47, color = "red")+
scale_alpha(guide="none")+
facet_wrap(~año, ncol=6)+
scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
guides(colour = guide_legend()) +
theme_bw()+
theme(panel.background = element_rect(fill = 'white'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
mas2
Falta componer mapas de tallas medias y alguna otra variable de interés.
CONCLUSION
En función de los datos analizados, el modelo propuesto es un enfoque de producción global estado espacio descrito por Pedersen & Berg (2017) y Acom (2015).
REFERENCIAS
Acom, I. C. M. (2015). ICES WKLIFE V REPORT 2015 Methodologies based on Life-history Traits , Development of Quantitative Assessment Exploitation Characteristics and other Report of the Fifth Workshop on the Relevant Parameters for Data-limited Stocks ( WKLIFE V ) Lisbon , Port. (October), 5–9.
Canales, C. M., Punt, A. E., & Mardones, M. (2021). Can a length-based pseudo-cohort analysis (LBPA) using multiple catch length-frequencies provide insight into population status in data-poor situations? Fisheries Research, 234(October 2020), 105810. https://doi.org/10.1016/j.fishres.2020.105810
Hordyk, A., Ono, K., Sainsbury, K., Loneragan, N., & Prince, J. (2014). Some explorations of the life history ratios to describe length composition, spawning-per-recruit, and the spawning potential ratio. ICES Journal of Marine Science, 72(1), 204–216. https://doi.org/10.1093/icesjms/fst235
Pedersen, M. W., & Berg, C. W. (2017). A stochastic surplus production model in continuous time. Fish and Fisheries, 18(2), 226–243. https://doi.org/10.1111/faf.12174
Rudd, M. B., & Thorson, J. T. (2018). Accounting for variable recruitment and fishing mortality in length-based stock assessments for data-limited fisheries. Canadian Journal of Fisheries and Aquatic Sciences, 75(7), 1019–1035. https://doi.org/10.1139/cjfas-2017-0143